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Upcoming gravitational wave-experiments promise a window for discovering new 
physics in astronomy. Detection sensitivity of the broadband laser interferomet- 
ric detectors LIGO/ VIRGO may be enhanced by matched filtering with accurate 
wave-form templates. Where analytic methods break down, we have to resort 
to numerical relativity, often in Hamiltonian or various hyperbolic formulations. 
Well-posed numerical relativity requires consistency with the elliptic constraints of 
energy and momentum conservation. We explore this using a choice of gauge in the 
future and a dynamical gauge in the past. Applied to a polarized Gowdy wave, this 
enables solving all ten vacuum Einstein equations. Evolution of the Schwarzschild 
metric in 3+1 and, more generally, sufficient conditions for well-posed numerical 
relativity continue to be open challenges. 



1 Introduction 

Laser Interferometric Gravitational- wave Observatories LIGO/ VIRGO 
is a broad band detector—targeting gravitational radiation from compact 
sources of a few solar massesE3. Notable sources are binary coalescence of neu- 
tron stars and black holes, as well as emissions from black hole-torus systems 
as recently proposedl£3. Their gravitational wave emissions provide a record 
of these strongly interacting catastrophic events and could contain most of the 
total energy released. What is the first source that LIGO/VIRGO may detect? 
Our knowledge of the gravitational wave-forms may be the determining factor 
in answering this question. For— this reason, numerical simulations of general 
relativity (numerical relativity e3) are receiving much attention in efforts to- 
wards matched filtering in searches for binary coalescence involving neutron 
stars and black holes Ell. 



1.1 Some astrophysical problems for numerical relativity 

There exists a broad spectrum of candidate sources of gravitational waves. 
Many sources have been proposed—for instance the coalescence of binaries of 
neutron stars .and-black holes c3e3£3, supernovae (seeEj), and rapidly spinning 
neutron stars&BEl. 

The gravitational wave-forms produced by neutron star-neutron star coa- 
lescence in the inspical phase are well-understood with post-Newtonian expan- 
sion technique (seelia for a recent review). Black hole-black hole coalescence 
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is very promising because of the expected larger amplitude signal. However, 
their event rate is highly uncertainty. Their wave-forms in the merger phase, 
believed to be relevant for LIGO/VIRGO detection, is noLpazEll understood. 
This is left as a challenge for numerical relativity (see, e.g.,ac3EJ). 

A new model was recenife. proposed for gravitational radiation from at 
torus around a black hole Ej'E3eZI. One particular feature is that it is asso- 
ciated with long gamma-ray bursts, whose intrinsic durations are about 20s 
on average. These energetic events are observed at a rate of about 2 per day 
and, corrected for beaming, at a rate of about 1 per year within a distance of 
lOOMpc. An important problem is the nature of the torus' mass-quadrupolc 
moment, which may be determined by numerical relativity. 

1.2 Numerical relativity: integration on thin ice 

Long-time integrations for the purpose of matched filtering requires accurate 
conservation of the elliptic constraints representing energy and momentum con- 
servation. This is a natural requirement, which becomes apparent in numerical 
experiments E3. 

The intrinsic hyperbolic structure of wave-motion permits well-posed ini- 
tial value problems, as follows from the energy method. In the continuum limit, 
conservation of energy and momentum is exact and drops out of the "energy 
balance sheet" . Consequently, well-posedness for general relativity reduces to 
the anticipated results for its hyperbolic structure. 

The nonlinear nature of the Einstein equations tends to introduce numer- 
ically a departure from exact conservation of the elliptic constraints. And the 
initial value problem for non-hyperbolic equations is often ill-posed, as for in- 
stance the Laplace and backward heat equation. The familiar Hadamard counter 
example to well-posedness is given by the solution 

u n = — sin nx sinh ny (1) 
of Laplace's equation on the upper half-plane: 

u X x + u yy = (— oo < x < oo, y > 0) (2) 
subject to the Cauchy data 

u(x,0)=0, u y (x, 0) = — sinncc. (3) 

The solution blows up in the face of large n, even though the initial data ap- 
proach zero in the norm of continuously differentiable functions. Numerically, 
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this tends to result in ill-poscdness: rapidly growing errors, regardless of the 
accuracy of initial conditions. The backward heat equation serves to illustrate 
similar ill-posedness in the presence of a first-order time-derivative,—. 

The observed instabilities associated with constraint violationsE3 suggests 
a link between variables with time-derivatives in the energy-momentum con- 
straints and ill-posedness. Well-posed numerical relativity requires these con- 
straints to be evolved in a consistent manner. A key test is the 'dynamical' 
evolution of a Schwarzschild black hole, where the dynamics derives from a 
singularity avoiding foliation of spacetime (see, e.g.,E3). In the absence of a 
covariant separation of the hyperbolic-elliptic structure of general relativity, 
we shall in this lecture discuss a recent proposal for an advanced-retarded evo- 
lution of the Einstein equations for exact preservation of the constraints under 
dynamical evolution. 

2 The Einstein equations, spacetime foliation and conservation laws 

The Einstein equations describe the structure of a curved spacetime manifold 
M with four-covariant metric g a b in response to a stress-energy tensor T ab . 
They are 

1 , . 

Rab - 2 gabR = 8%T ab (4) 

as an equation for the Ricci tensor = R c bcd and its scalar curvature R = 
R c c . These are expressions in terms of the Riemann tensor R a bcd - The left 
hand-side is commonly denoted by the Einstein tensor G ab , with the property 
that V°G a b = (the Bianchi identity) is consistent with energy-momentum 
conservation V a T af> = 0. We are at liberty to include a cosmological constant 
term —Ag a i, on the right hand-side of (0). The equations of motion (|J) derive 
from the Hilbert action 

S = ( y^Rd^x. (5) 

J M 

Translation invariance is a symmetry in this action (||). By Noether's theorem, 
this leads to four conservation laws of energy and momentum. The associated 
gauge group is the Lorentz group SO(3,l,i?) of boosts and rotations. These 
conservation laws become explicit on spacelike hypersurfaces Y, t : t =const., 
where t denotes a timelike coordinate. Combined, the surfaces T, t provide a 
foliation of spacetime. 
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2. 1 Spacetime foliation in spacelike hypersurfaces E t 
Hypcrsurfaces E t of constant time come with two vectors: 



Na=gta, n a = d a t/^-d a td a t, (6) 

where n a denotes the_unit normal (n 2 = —1) to E t . (The vector J\f a is com- 
monly denoted by i a E3.) Generally, the covariant vectors Af a and n a are inde- 
pendent. Marching from one hypersurface to the next brings along a variation 
dt, along with the covariant displacement 

ds a = N a dt. (7) 

The displacement ds^ expresses M a as a "flow of time." It can be expressed 
in terms of orthogonal projections on n a and S t , in terms of the lapse function 
N and shift functions N a , 

M a = Nn a + N a . (8) 

Here TV — —AT a n a and N a = h^Afb, expressed in the metric 

hab = gab + n a n b (9) 

as the orthogonal projection of g a b onto Et. Note that ds 2 = JV 2 dt 2 = gttdt 2 
as the square of (0), so that g tt = —N 2 + N C N C . With n a = (n t , 0,0,0), it 
follows that 

(N C N C -N 2 NA 

Nl (10) 

where i,j refer to the spatial coordinates x l of (t,x l ). The lapse function 
satisfies \J—g = NVh. The four degrees of freedom in the five functions 
(N,N a ) are algebraically equivalen±_to Af a . An equivalent expression for the 
line-element, in so-called 3+1 formE3, is 

ds 2 = -a 2 dt 2 + j ij (dx i + f3 l )(dx j + /3 j ), (11) 

where a = N is referred to as the redshift factor and "fij/3 J = gu- 



2.2 Conservation of energy and momentum 

Coordinate invariance introduces a certain degeneracy in the Einstein equa- 
tions. There are no second-order time-derivatives of g a b in the components 
G n b = G a bn a of the Einstein tensor G a b — Rab — (l/2)g a bR- Consequently, 
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the expression G n b forms entirely out of Cauchy data on S t (data and their 
first time-derivatives). The embedding of £t in four-dimensional spacetime 
is expressed in terms of the symmetric extrinsic curvature tensor K a \,. If fib 
denotes a unit tangent to a geodesic orthogonal to S t , then 

K a b = V a n b = ^C n h a b- (12) 

Thus, K a b represents a time-like derivative of the—metric in £ t , which is a 
velocity of h a b- We use here the sign convention inEa; K a b with opposite sign 
is commonly used in numerical relativity. We then have 

D b K b a - D a K = 8nT an , 
WR + K 2 -K ab K ab = l6irT nn [U) 

a consequence of|the-QK>jection ^R a bcd of the four-dimensional Riemann tensor 
R a bcd on E t (seeO'ElEil for detailed calculations). Here, K = K£ denotes the 
trace of K. These expressions ( |l3| ) are, respectively, the conservation laws 
of linear momentum and energy. These equations are elliptic in the spatial 
coordinates internal to £(. 



3 Two marching methods for hyperbolic formulations 

A practical frame- work for numerical relativity consists of marching data from 
one space-like hypersurface S t : t =const. to the next. The hypersurface E f is 
generally dynamical. A slicing of spacetime either comes before or after a choice 
of dynamical variables. Chosen before, the variables live in £4 as projections 
of the underlying four-covariant metric. Chosen after, one continues to work 
with four-covariant metric parametrized over the hypersurfaces T, t - There is 
no dictum for the order of these choices, but they do give manifestly different 
formulations. 



3. 1 Slice first: the Hamiltonian approach 

In the Hamiltonian approach, we consider first a choice of foliation of spacetime 
in spacelike hypersurface S t of constant coordinate time t. Their dynamics is 
described by the projected metric h a b with canonical momentum TT a b, satisfying 
the Hamiltonian equations 
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where the dot denotes the Lie derivative £ t with respect to the vector field 
t a = M a of the flow of time (|J). This Lie derivative reduces to differentiation 
with respect to t in our coordinate system {t,x l ). 

The Hamiltonian equations derive from the Hilbert action (||) . An excellent 
presentation is given in Appendix E inl£3, which is briefly summarized here. 

The Lagrangian density C = R^—g can be expressed as a sum of the 
three-curvature of S t and a quadratic form of the extrinsic curvature 
tensor K ab = V a n 6 , given by C = VhN[^R + K ab K ab - K 2 }. We haveS 
K ab = L n h ab /2 = [h ab - D a N b - D b N a ]/(2N), where D a = h b a V a denotes the 
derivative internal to £*. It follows that 

i: ^ = '±-=yfh(K ab -Kh ab ). (15) 
dh 

Coordinate invariance of the Einstein equations leaves gta and, hence, 
(N,N a ) freely specifyable. Tracing back, we indeed find no first-order time 
derivatives of gtb in the Lagrangian density £, whereby the associated canoni- 
cal momenta vanish. The Hamiltonian density associated with the dynamical 
degrees of freedom reads therefore H = ir ab h ab — C. The variational derivative 
of H — J s TL\fhd i x with respect to these gauge functions obtains the con- 
servation laws of energy and momentum (E~3[). The variational derivative with 
respect to the dynamical variables (h ab , 7r a ^obtains the ADM formulationEEj 

h ab - 2D (a N b) = 2h- 1 / 2 N(TT ab - h ab n) (16) 

and 

n ab - 2ir< a D c N b '> = -Nh^ 2 {^R- ± (3) ' Rh ab ) 

+h 1 ' 2 {D a D b N - h ab D c D c N) 

+ \Nh- 1 / 2 h ab {^:^-\^ 2 ) (17) 

-2A^/l- 1 / 2 (7T aC 7rJ - ±7T7T afc ) 

+h 1 / 2 D c (h- 1 / 2 N c n ab ). 

In numerical relativity, these Hamiltonian evolution equations are often 
considered in terms of the non-canonical pair (7^- , Kij ) , with 7^ as in ([llj) and 
Kij, where i,j refer to the spatial components in (t,x % ). Thus, ( |l6| ) and (JTt[) 
become, using the vacuum case of (|l3|), 



(dt - Cp)jij = ~2aKij, 

(d t - £p)K tJ =-D l D J a + a[^R l3 - 2Kf j + KKy], 



(18) 



where we use the definition of Kij with opposite sign of (|12|), following the 
convention in this context. 
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3.2 Hyperbolic formulations in the Hamiltonian approach 



Hyperbolic systems of equations in the non-canonical variables (7^ , Kij ) , or 



closely if^at 
notably! 3 



aye been derived in various forms by several groups, 
seeEU for a comprehensive review. This approach 
typically comprises constraints on the lapse functions. For a recent comparison 
study between fopoyi^ations in Hamiltonian variables and related hyperbolic 
Different formulations display various degrees of nu- 



formulations, see 
merical stability!^ ^ 
of hyperbolicityEJO 



is 



which appears not to depend significantly on the degree 



3.3 Dynamical conservation of constraints 

The evolution equations for general relativity may formally be modified, such 
that the energy and momentum constraints become a stable manifold of phys= 
ical solutions. This has recently: been considered in a linearized treatments 
and in Ashtekar's formulationEjOo. This holds some promise in providing a 
unified treatment of the dynamical and the elliptic parts of general relativity. 
Numerical results on accuracy and stability are inconclusive at present E3. 



3.4 Slice last: the four-covariant approach 

General relativity can be written as nonlinear wave equations for the Ricmarm- 
Cartan connections in the tetrad formulation. This builds on Pirani's-argu- 
ments concerning the role of the Riemann tensor in gravitational wanes E3 and 
on Yang-Mills formulations of generaLxel ati vity, following UtiyamaEa and de- 
veloped by Ashtekar and co-workers EJOEll. Starting point is a divergence 
equation for the Riemann tensor with an anti-symmetric derivative of the 
stress-energy tensor as a source-term. 

The interwoveness of wave motion and causal structure distinguishes grav- 
ity from other field theories. This becomes apparent in nonlinear wave equa- 
tions for the connections on the curved spacetime manifold side-by-side with 
equations of structure for th e e volution of the metric in the tangent bundle. 

The tetrad approach EJE3 bears some relation to but is different from 
Ashtekar's propram on nonperturbative quantum gravity. The original Ashtekar 
variables are SU(2,C) soldering forms and an associated complex connection in 
which the constraint equations become polynomial. The Riemann-Cartan vari- 
able is a real SO(3,l,R) connection. In Ashtekar's variables, a real spacetime is 
recovered from the complex one by reality constraints. See Barbero for a trans- 
la-tion of Ashtekar's approachinto SO(3,R) phase space with real connections 
I13il3. The main innovation infill is the incorporation of the Lorentz gauge con- 
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dition ( Pq ) which obtains new hyperbolic evolution equations in four-covariant 
form (below). 

Following Pirani, we take the view that gravitational wave-motion is con- 
tained in the Riemann tensor, R a bcd- It satisfies the Bianchi identity 

3V[ e i? Q b] c rf = VeRabcd + V a Rbacd + V ' bReacd = 0. (19) 

This gives rise to the homogeneous divergence equation V a * R a bcd — 0, where 
*Rabcd = (1/2)Cq6 Refcd denotes its dual. Upon interaction with matter in 
accord with the Einstein equations, the Ricci tensor satisfies R a b — 8ir[T a i, — 
\g a bT\. The Bianchi identity above also gives V d R a bcd = 2V[ h i? a ] c , which 
obtains the inhomogeneous divergence equation 

V a Rabcd = l6n(V [c T d]b - ~<? 6[d V c] T). (20) 

The quantity on the right hand-side shall be referred to as \QirT bcc i. This term is 
divergence free, V b Tb C( j = 0, on account of the conservation law V a T ab = and 
consistent with divergence-free condition V b V a R a bcd = on the left hand-side 
( ^p| ) (by anti-symmetry of the Riemann tensor in its first two indices). 

Introduce the Riemann-Cartan connections u) aflI/ = (e M ) c V a (e u ) c associ- 
ated with a tetrad {(e p ) 6 }^ =1 . Then the above-mentioned homogeneous and 
inhomogeneous divergence equations take the form 

V a * R ab ^ = 0, V Q i? a &^ = I6n^, (21) 

where the w ^„ define a gauge-covariant derivative in accord with the Yang- 
Mills construction V a = V a + [ui a , ■}. The first of ( f2l| ) gives rise to the represen- 
tation Rab^v = V a a;b MI , - Vb^a/Mu + [w a ,Ub]n V . The gauge covariant derivative 
satisfies the identity V a (e A1 ) b = 0, which implies the equations of structure 
<9[ a (e P )b]) = (e")[i. w a]^ " leaving d t (e M )t undefined. Next, define £ b = (<9 t ) b , 
and introduce the tetrad lapse functions 

N» = (e^) a C (22) 

as freely specifiable functions. Thus, the equations of structure become a 
system of first-order ordinary differential equations 

d t (e^) b + <(e„) b = dbN M + uj b »N v . (23) 

The tetrad lapse functions are algebraically equivalent to the familiar Hamil- 
tonian lapse, N, and shift functions, N a , through (|l0|): 

g a t = N a (e a ) a = {N q N q - N 2 ,N p ). (24) 
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The term ujb llu N u on the right hand-side of ( |23| ) shows that the tetrad lapse 
functions introduce different transformations on each of the legs; the term 
uj' t ^{e v )b on the left hand-side introduces a transformation with applies to all 
four legs simultaneously. It follows that it is the infinitesimal Lorentz trans- 
formations uj' t which provide the internal gauge transformations. 



3.5 Hyperbolic equations in the four-covariant approach 

An important aspect here is internal gauge fixing on the Lorentz group SO(3,l,i?) 
associated with the choice of tetrad. To fix gauge, we propose using the Lorentz 
gauge Ed 

Cfiu := V a uj atlu = 0. (25) 

This fixes unique evolution equations for the internal gauge, resulting in nonlin- 
ear wave equations for the connections u) aiiV . These complement the equations 
of structure for the evolution of the tetrad legs, and together form a complete 
system of evolution equations. The Lorentz gauge (|5|) defines a choice of ac- 
celeration of the tetrad legs, through the infinitesimal Lorentz transformations 
ujt^v mentioned above. In a different context of compact gauge groups and 
a metric with Euclidean signature, its-geometric significance has been inter- 
preted by Lewandowski et al. (1983) cil. These six constraints (p5|) can be 
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a hyperbolic implementation by application of the divergence technique 



^""{Rab^u + QabC^u} = 167TT fcpIy , (26) 



which preserves c^ v = is preserved in the future domain of dependence of 
the support of physical initial data. By explicit calculation, we have 

Qo; 0jUI / - RaU cf j.v - [^ C , V Wc]/w = 16-KTafj.u, (27) 

where □ denotes the Yang-Mills wave operator V 2 . The Ricci tensor on the 
left hand-side may be understood in terms of T a b using_thc Einstein equations. 
The above provides the following covariant separationEiJ. 



Theorem 1. Gravitational waves propagate on a curved spacetime manifold 
by nonlinear wave equations in a Lorentz gauge on the Riemann-Cartan con- 
nections. In response to this wave motion, the causal structure of the manifold 
evolves in the tangent bundle by the equations of structure. The Hamiltonian 
lapse and shift functions find their algebraic counterparts in the tetrad lapse 
functions N^. 
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We remark that away from the matter source, the vacuum equations read 

Ow a/n ,-[w c ,V a WcW = 0. (28) 

This .siacuum case has been considered numerically in one-dimensional Gowdy- 
testsE3 by using the underlying first-order system for the components of the 
Riemann tensor (p6|). 



4 Hyperbolic Einstein-MHD equations 

Relativistic hydrodynamics and magnetohydrodynamics has received consid- 
erable attention in the simulations of astrophysical jetscSEH. Recentbz. these 
efforts result in simulations of astrophysical jets around black holes E3 with 

a tensions to flows around rotating black holes for a few dynamical time-scales 
. The latter shows a transition of accretion disk outflows towards a state of 
differential rotation in the vicinity of the black hole. 

GRBs from rotating black holes are associated with a_compact torus or 
disk, representing-hinary black hole-neutron star coalescenceo, failed-supernovae 
EH or hypernovaeocj. The torus or disk may well be magnetized with the rem- 
nant field of the progenitor star, i.e., the neutron star in the coalescence sce- 
nario or young massive star in the failed-supernova or hypernova scenario. The 
suggests simulating the creation of gravitational waves by high-density matter 
around a black hole in the approximation of ideal magnetohydrodynamics. 

A hyperbolic formulation of ideal magnetohydrodynamics, used in the sim- 
ulation of relativistic jetsLil, is 

V a T ab = 0, 

V a (uI Q /i b l + g ab c) = 0, (29) 
V a {ru a ) = 0, 

expressing conservation of energy-momentum, Faraday's-equations in diver- 
gence form which conserves the constraint c = u c h c = EHI, and conservation 
of baryon number. This ia,a single fluid description of an ideal, inviscid fluid 
with stress-energy tensorca 

T ab = (r + 7^7 (7 - 1) + h 2 )u a u h + (P + h 2 /2)g ab - h a h\ (30) 

where u a denotes the velocity four-vector, r the comoving rest mass density 
and P the pressure with polytropic equation of state P = Kr 1 and polytropic 
index 7. 

Perhaps Theorem 1 and ( p^ ) may serve as a starting point for hyperbolic 
Einstein-MHD equations for the purpose of numerical simulations. 
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5 Past and future gauge in numerical relativity 

The Einstein equations pose six equations for dynamical evolution plus the 
four constraints of energy and momentum conservation. The latter involve 
the normal vector n a to the surfaces of foliation, i.e.: the projection operator 
onto these surfaces. The six dynamical degrees of freedom pertain to variables 
subject to second-order time-derivatives, while gauge-variables are subject to 
constraints on their first-order time-derivatives. In a discrete setting, the first 
live on three and the second on two hypersurfaces of constant time. This 
suggests to consider the ten degrees of freedom involved to be the six dynam- 
ical degrees of freedom supplemented with four gauge-variables from the past. 
The gain is exact conservation of energy-momentum, traded off against exact 
projections in the past. 

Non-exact projections naturally permit an uncertainty between the three- 
metric and its canonical momentum within the underlying context of a four- 
covariant theory, i.e.: also in regards ±o the association with the hypersurfacc 
at hand. In the covariant approach offij, this would thus reflect an uncertainty 
in the tetrad elements, which define the projection, and their connections. This 
points towards a potential connection to quantum gravity. Indeed, soon after 
this work was proposedEa, the author learned of a very interesting independent 
discussion on the problem of consistent discretizations in this contextEZl. 

5.1 A discretized initial value problem 

We illustrate our this approach on the vacuum Einstein equations 

Rab = 0. (31) 

The Ricci tensor R ao is a second-order expression in the metric g a b, whereby 
( |3l| ) defines a relationship between metric data , g" b , g^ 1 ) on a triple of 

time-slices t n —\ < t n < t n+ \: 

Rab(g n ab + \g: b ,9 , a l b 1 )=0- (32) 

Here, Rm = R a bc A derived from the Riemann tensor 

R a bed — dd^bc ~ 9 c ^bd + ^bc^ed ~ ^bdFec- (33) 

This expression ( [33] ) can be discretized by finite differencing on a triple of 
time-slices with preservation of the quasi-linear second-order structure of R a b- 
Algebraic gauge-fixing takes the form of specifying the components iV a = 
gta in coordinates {x a }* =1 with t = x 1 time-like. A gauge-choice on a triple 
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of time-slices amounts to a choice of (iV™ -1 , N™, iV™ +1 ). This gauge-choice in 
the metric arises explicitly in the Hamiltonian constraints of energy-momentum 
conservation. The components hij — gij , where i,j = 2,3,4 refer to projections 
of the metric into the time-slice t =const., which describe the dynamical part 
of the metric. The combination (hq, N a ) reflects the mixed hyperbolic-elliptic 
structure in numerical relativity and ( |3l|) represents ten evolution equations 
in these variables on a triple of time-slices. 

In algebraic gauge- fixing, we prescribe -/V™ +1 as a future gauge in com- 
puting on a future hypersurface t = t n +\ from data at present and past 
hypersurfaces t — t n _i and t — t n . Re-introducing -/V™ -1 as dynamical gauge 
in the past gives closure, leaving h™~ x fixed. This combination of ten degrees 
of freedom defines an advanced hyperbolic-retarded elliptic evolution of the 
metric. The paritioning of the metric in past and future variables as 



The dots refer to the remaining data (ft.™- , hy, -/V™, JV™ +1 ), which are kept 
fixed while solving for AT™ -1 ). Thus, ( |35| ) which takes into account 

all ten Einstein equations with no reduction of variables. Time-stepping by 
( |35| ) evolves the metric into the future with dynamical gauge in the past, in 
an effort to satisfy energy-momentum conservation within the definition of the 
discretized Ricci tensor. Because (^) comprises derivatives of N a only to first- 
order in time, numerically though the data N™ +1 and N™~ 1 , we anticipate that 
the evolution of N a is of first-order in the t— discretization At. This introduces 
non-exactness in h™j as projections of g a t, on t — At to within the same 
order of accuracy. It may result in a first-order drift in the t— labeling of the 
hypersurfaces - permitted by coordinate invariance. 

5.2 A polarized Gowdy wave 

The presented approach can be illustrated on a polarized Gowdy wave. Gowdy 
cosmologies have compact space-like hypersurfaces with two Killing vectors <9 CT 
and ds ■ With cyclic boundary conditions, the space-like hypersurfaces are 




(34) 



thus obtains ten dynamical variables in the ten equations 



R ab (hff\NZ-\---)=0 at t = t 



(35) 
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homcomorphic to the three-torus as a model universe collapsing into a singu- 
larity. The associated line-element is (see, e.g. ,113) 

ds 2 = e^-^' 2 (-e- 2T dr 2 + d9 2 ) + d£ 2 , (36) 

where A = A(t, 0) and dH denotes the surface element in the space supported 
by the two Killing vectors. Polarized Gowdy waves form a special case, which 
permit a reduction to 

d£ 2 = e- T (e p da 2 + e - p dS 2 ) . (37) 

Here P satisfies a linear wave-equation P TT = e~ 2T Pee; a long wave-length 
solution is 

Po(T,0) = Y o (e- T )cos6, (38) 
where Yq is the Bessel function of the second kind of order zero. This leaves 

AM) = iy o (e- T )yi(e- r ) e - r cos20 + i^ ^ (F ' 2 ( S ) + Y 2 (s)) sds. (39) 

A spectrally accurate numerical integration is described in El. 

The implicit equation ( |35| ) for the dynamical variables (h™^ 1 , N™^ 1 ) has 
been implemented numerically. We have done so by solving for the all ten 
components (/i™ +1 , N^ 1 ) using Newton iterations on these variables. This 
procedure uses a numerical evaluation of the Jacobian 

where the capital indices A,B = 1, 2, • • • , 10 refer to the labeling 

Ra = (Rn, R22, R33, R44, R12, R13, R14, -R23, -R24, ^34), / 1 1 \ 



The Ricci tensor (|33|) has been implemented by second-order finite differencing, 
such that it remains quasi- linear in the second derivatives. In particular, the 
Christoffel symbols 

Kb = \g Ce (9cb,a + 9ac,b - 9ab,c) (42) 

is obtained by symmetric finite-differencing on the metric components, and 
itself differentiated by the product rule following individual numerical differ- 
entiations of g ab and (g c b,a + g ac ,b — g a b,c)- The choice of future gauge iV™ +1 is 
provided by the the components 

g at = (e^-W, 0,0.0) (43) 
13 



of the analytical line-element (36-39|), which facilitates error analysis by direct 
comparison of the numerical results with the analytic solution. It will be 
appreciated that in principle other choices of N£ +1 can be made. 

Fig. 1 shows numerical results for evolution of initial data on the interval 
< r < 4. The results show that all Einstein equations in the form of R a b = 
are satisfied with arbitrary precision, while the metric components are solved 
accurately to within one percent. The asymptotic behavior of the implicit cor- 
rections to the lapse functions are shown in Fig. 2. Note that these corrections 
are finite to first-order in At: the corrections SN on the past lapse satisfy 

= Nn+1 ^)- Nn - 1 ^) = o (Ar) (44) 

This first-order dependence is a testimonial to the fact that the lapse function 
appears in the Einstein equations to first-order in time. 



Theorem 2. A choice of gauge in the future and a dynamical gauge in the 
past obtains a discretization of general relativity consistent with energy and 
momentum conservation, permitting all ten Einstein equations to be solved in 
ten dynamical variables in the presence on non-exact projections of the four- 
covariant metric on the past hypersurface of constant time. 



6 Summary and conclusion 

Well-posed numerical relativity is a long-standing challenge in the calculation 
of wave-forms for astrophysical sources of gravitational radiation. A necessary 
condition for stable numerical relativity is accurate conservation of the energy 
and momentum constraints ( "integration on thin W") .-^This has been pursued 
by implementing these constraints dynamically EjO'Ej Here, we have ex- 
plored a consistent discretization for exact conservation of energy-momentum 
constraints using a choice of gauge in the future and a dynamical gauge in the 
past. This permits integration of all ten Einstein equations, while allowing for 
in-exact projections of the four-covariant metric onto the surfaces of foliation 
of spacetime. The simulation of a nonlinear one-dimensional Gowdy wave by 
implicit time-stepping according to the ten discretized vacuum Einstein equa- 
tions (^) serves to illustrate a numerical implementation. 

A major open problem is obtaining sufficient conditions for stability. In 
this presented approach, it becomes of interest to consider novel definitions of 
the future gauge as a function of present gauge. We leave this as a suggested 
direction for future-development. In light of the recent discussion by Gambini 
and Pullin (2002) E3, the question arises: is well-posed numerical relativity 
related to consistent discretization in quantum gravity? 
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Figure 1: Shown is the simulation for < r < 4 of the polarized Gowdy wave. The 
solutions P(t, 8) and A(r, 8) are displayed as a function of (r, 8) (upper windows). The middle 
windows display the solutions for r = 4, wherein the circles denote the numerical solution 
and the solid lines the analytical solution. The r— evolution of the errors (lower windows) 
are computed relative to the analytical solution to Gowdy's reduced wave equation. The 
simulations use a discretization of 8 by mi = 64 points and the r— interval by m,2 = 1024 
time-steps. Particular to the proposed numerical algorithm is a dynamical gauge in the past 
and a prescribed gauge in the future. This permits satisfying all of the discretized Einstein 
equations R a i, = to within arbitrary precision by Newton iterations. The slight increase 
in the error of about 10~ 10 reflects the exponential growth of the analytic solution, because 
the Gowdy cosmology evolves towards a singularity. (Reprinted from van Putten, M.H.P.M., 
Class. Quantum Grav., 19, L51, ©IOP 2002.) 
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Figure 2: Shown are the self-consistent corrections on the slicing t = t n +l, introduced 
by the difference between the past gauge JV n_1 (t„_|_2) to the hypersurfacc t = t n +2 and 
the earlier future gauge N n+1 (t n ) to the hypersurfacc t = t n . The three curves refer to 
different discretizations mi = 16, 32 and 64 points with, respectively, m,2 = 256, 512 and 1024 
time-steps. These similar results for various discretizations indicate asymptotic behavior 
consistent with the first-order appearance of the lapse function in the Einstein equations. 
A first-order accuracy in lapse introduces a commensurate offset in slicing or, cquivalcntly, 
an offset in the coordinate t. Similar results obtain for the same spatial discretizations 
mi = 16,32 and m\ = 64 with time-steps at one-half the value, i.e., using r«2 = 512, 1024 
and, respectively, m,2 = 2048 time-steps. (Reprinted from van Putten, M.H.P.M., Class. 
Quantum Grav., 19, L51, ©IOP 2002.) 
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Evolving eternal Schwarzscl. 
problem for these developments t 



.d black holes in 3+1 may serve as a test 
1 More generally, it would be of interest to 
consider exact conservation of energy-momentum in higher dimensions, per- 
haps using a combination of any of the modern hyperbolic formulations and 
efficient elliptic solvers. 
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